/**
 * @file VC.h
 * @brief Implementation of V-cycle.
 * @author XDDDD
 * @version 
 * @date 2021-06-10
 */

#ifndef __PAISLEYPARK__VC_H__
#define __PAISLEYPARK__VC_H__

#include "MultigridSolver.h"

template <size_t DIM>
class VC : public MultigridSolver<DIM>{
public:
	VC() : MultigridSolver<DIM>(){};
	~VC(){};
	VC(std::string _RestrictionId, std::string _InerpolationId) : MultigridSolver<DIM>(_RestrictionId, _InerpolationId){};
	void Solve(VectorFunction<DIM>& _f, VectorFunction<DIM>& _g, std::vector<double>& _v,
			  const size_t _N, const size_t _M, const size_t itr_num, const size_t nu1, const size_t nu2, const double _epsilon) override;
};

template <size_t DIM>
void VC<DIM>::Solve(VectorFunction<DIM>& _f, VectorFunction<DIM>& _g, std::vector<double>& _v,
					const size_t _N, const size_t _M, const size_t itr_num, const size_t _nu1, const size_t _nu2, const double _epsilon) {
	bool is_stop = false;
	this->set(_v, _N, _f, _g);
	size_t N = this->__U.get_N();
	size_t pre_N = N;
	for(size_t j = 0; j != itr_num; j++) {
		is_stop = false;
		while(!is_stop) {
			for(size_t i = 0; i != _nu1; i++) {
					this->smooth(pre_N);
			}
			if(pre_N == _M) {
				is_stop = true;
				break;
			}
			this->set_e(pre_N);
			
			pre_N -= 1;
			std::vector<size_t> coord(DIM, 1);
			size_t pole = 0;
			while(pole != DIM) {
				pole = 0;
				if(fabs(this->get_F(coord, pre_N)) < _epsilon) {
					is_stop = true;
					break;
				}
				coord[pole]++;
				while(coord[pole] = (1 << pre_N)) {
					coord[pole] = 1;
					pole++;
					if(pole == DIM)
						break;
					coord[pole]++;
				}
			}
			this->__U.clean(pre_N);
		}
		if(pre_N == _M)
			this->solve_coarest(pre_N);
		else {
			for(size_t i = 0; i != _nu2; i++) {
				this->smooth(pre_N);
			}
		}
		is_stop = false;
		while(!is_stop) {
			this->putback_e(pre_N);
			pre_N++;
			for(size_t i = 0; i != _nu2; i++) {
				this->smooth(pre_N);
			}
			if(pre_N == N) {
				is_stop = true;
			}
		}
	}
};
#else
//Do nothing.
#endif
